Classification of early age facial growth pattern and identification of the genetic basis in two Korean populations

Childhood to adolescence is an accelerated growth period, and genetic features can influence differences of individual growth patterns. In this study, we examined the genetic basis of early age facial growth (EAFG) patterns. Facial shape phenotypes were defined using facial landmark distances, identifying five growth patterns: continued-decrease, decrease-to-increase, constant, increase-to-decrease, and continued-increase. We conducted genome-wide association studies (GWAS) for 10 horizontal and 11 vertical phenotypes. The most significant association for horizontal phenotypes was rs610831 (TRIM29; β = 0.92, p-value = 1.9 × 10−9) and for vertical phenotypes was rs6898746 (ZSWIM6; β = 0.1103, p-value = 2.5 × 10−8). It is highly correlated with genes already reported for facial growth. This study is the first to classify and characterize facial growth patterns and related genetic polymorphisms.

Defining each measurement and analyzing the time series of facial measurements according to age. The measurements from each past photograph were compared against those of the current photograph using non-linear regression methods to determine the facial changes over time. The facial growth patterns of each individual are determined by the visual inspection of individual's changes were used for the genetic association study.
To understand the changes in facial measurements with age, a graph of the time series for each individual facial measurement was plotted for 21 phenotypes according to the age relative to the current age using a nonlinear model. Supplementary Fig. S1A-C graphically represent individual growth patterns for each representative eye, nose, and mouth phenotype. Table 1 summarizes the distribution of facial growth patterns by face region in Pop1 and Pop2. We identified five EAFG patterns: Pattern 1 (DD), continued decrease; Pattern 2 (DI), decrease to increase; Pattern 3 (CC), constant; Pattern 4 (ID), increase to decrease; and Pattern 5 (II), continued increase (Fig. 2). The X-axis represents age, and the Y-axis represents facial distances between two selected points www.nature.com/scientificreports/ based on the 19 features measured. Among the EAFG patterns, the higher the frequency, the darker the gray. Among 21 phenotypes, 14 phenotypes, H1-H7 and V5-V11 showed similar distribution patterns between POP1 and POP2, whereas 7 phenotypes showed unique distribution patterns between POP1 and POP2. However, the 7 phenotypes representing unique growth patterns were clustered within Pattern 1 and Pattern 5. High proportions of Pattern 5 were observed for both horizontal and vertical measurements. When examining specific regions of the face, the area around the eyes showed high proportions of Pattern 1, whereas the other facial phenotypes showed high proportions of Pattern 5. For most Koreans, the distance between H1 and H2 decreased with aging, the measurement for H3 increased with age, and H4 showed a tendency to decrease with age. Among the other eye widths, H5 and H6, most commonly increased with age, and the vertical eye measurements, V1 and V2, also showed a tendency toward an increasing pattern. In the nose, the width of the nostrils most commonly increased, and the length of the vertical axis of the nose showed the largest increasing pattern. Around the lips, many individual's patterns showed measurements that were maintained or increased types as they aging.
Genotype analysis. Genome-wide single-nucleotide polymorphism (SNP) genotypes were obtained from an 800 K SNP microarray experiment using an Axiom array followed by imputation using 1000 Genomes Phase 3 data 18 , resulting in a total of 7,375,270 polymorphic SNPs included in the GWAS. We conducted a GWAS for the combined analysis of POP1 and POP2, in addition to separate analyses for POP1 and POP2. The significant or suggestive SNPs from the combined analysis were determined based on the criteria of a p-value < 5 × 10 −8 for significant SNPs and 5 × 10 −8 ≤ p-value < 1 × 10 −5 for suggestive SNPs. In addition, SNPs in the individual POP1 and POP2 analysis with p-value < 0.05 were considered significant. The combined GWAS results are illustrated using quantile-quantile (QQ) plots ( Supplementary Fig. S2) and Manhattan plots ( Supplementary Fig. S3) for each phenotype. A total of 97 SNPs satisfied the genome-wide significance criteria (p-values < 5 × 10 −8 ), and 759 SNPs were identified with suggestive association p-values (5 × 10 −8 ≤ p-value < 1 × 10 −5 ) for 21 facial phenotypes (Supplementary Table S3). The SNPs were analyzed by clustering patterns: 77 SNPs were singletons (i.e., a single significant SNP without co-segregation with other SNPs nearby the significant SNP), and 729 SNPs in 104 loci showed a clustered pattern (i.e., significant or suggestive SNPs co-segregated with more than three other significant or suggestive SNPs). The significant or suggestive and clustered SNP loci were illustrated using regional signal plots ( Supplementary Fig. S4). The top significant SNPs for each significant cluster are described in Tables 2 and 3. The criteria for dividing DD, DI, CC, ID, and II patterns were established and validated using quantitative trait association analysis (Wald test). Among the 5 EAFG patterns in the horizontal and the vertical measurements, we identified significant and suggestive SNPs. The phenotype used in this study was analyzed by coding facial growth patterns from the past to the present from 1 to 5 (Supplementary Table S4). There was no significant difference in the distribution of patterns between males and females in the facial growth pattern. Therefore, gender was not used as a covariate in the GWAS analysis in this study. The SNPs of the existing Facial Measurement GWAS results were checked once more in the results of this study, and the results were included in a separate Supplementary Table S5. We could test 20 SNPs that was  www.nature.com/scientificreports/ previously reported in other studies, and described the association results for the early facial growth patterns that the major phenotypes of this study. Most of the tested SNPs were found to have no significance or very weak significance. The reason for this difference is that this study is not an analysis to discover genetic markers related to differences in facial shape between individuals, but an analysis of childhood facial growth patterns. A limitation of this study was that there were not enough participants and some factors that could affect facial growth were not measured. The facial measurements were appeared 2D images rather than 3D images. To understand for facial growth patterns, we tried to measure for facial traits, such as take a photograph of early age participants in the studio and recruited past pictures. Through like this study, further study will be more understand for facial morphology and facial growth.

Discussion
Human development is characterized by distinct developmental processes, especially during adolescence, and the speed and direction of craniofacial development differ for each person. Therefore, under the assumption that the speed and direction of development would differ across individuals, we calculated the changes in facial measurements according to age over time for each individual and for each indicator and divided these change patterns into five major categories. An index was calculated based on the positioning of landmarks in the facial profile picture, according to a widely used method for current facial analyses. GWAS analysis of facial development patterns was performed by recoding each of the five growth patterns as individual values.
Despite the significant difference from the existing approach such as GWAS analysis of face measurements or facial deformities, we were able to obtain GWAS results that were repeatedly associated with facial features. A total of 97 significant indicators were identified, including indicators related to craniofacial development in 19 areas. Because the probabilistic effects and differences of facial phenotypes must be confirmed by replication analysis of identified loci, we divided the collection into two groups, and the genetic influences in each group were analyzed for stochastic effects through replication analysis and it were confirmed through the replication indicators.
In the nose, the width of the nostrils most commonly increased with age, and the length of the nose showed the largest tendency toward increase across all of the vertical axis. Around the lips, many individuals showed patterns of increase or maintenance, whereas the patterns associated with the mouth were distributed differently in each population, indicating the degree of difference in the growth patterns among individuals.
The length of H1 and H2 have the most decreasing patterns with aging. H5 and H6 have the most increasing patterns with aging. The most noticeable changes were observed for the cutis and subcutaneous bone.
For both POP1 and POP2, as shown in Table 1, H3 most commonly have the most increasing patterns, whereas H4 most commonly have the most decreasing patterns. Because as they grow, the elasticity of adjacent tissues under the eyes decreases due to growth, resulting in the sagging tail appearance of the eyes. In addition, H8, H9, and H10 appeared to be very significant indicators, which appeared to influence each other. H8 is a  www.nature.com/scientificreports/ diagonal length on the left side of the nose, which tends to increase with age. H9 is the width of the nose, which tends to increase with age because the lower lateral cartilage and the skin surrounding the ends of the septum weaken, losing elasticity. Among the vertical axis lengths, many indicators showed similar patterns of increase as the horizontal area, and the vertical axis indicators of the face appear to affect each other during growth. The vertical lengths of the eyes, V1 and V2, showed the greatest tendency toward an increasing pattern. The length of the nose, measured by V5-V8, commonly increases until the age of 20 years. The characteristics of the nose are well known and include major changes, such as long, drooping tips 19 . The bone base that supports the nose in youth, a pair of  www.nature.com/scientificreports/ nasal bones, and the ascending process of the maxilla are responsible for many of the soft tissue changes that are observed in the nose during aging 20 . The pattern frequencies measured for EAFG showed that although we used independent populations, our results were replicated in each population. As shown in Table 1, approximately 70% of the facial development patterns were replicated in each group. The index with the highest frequency was replicated in each group, indicating a common pattern across populations. Although a few indices showed a different pattern, these unique indices clustered into large categories (increase or decrease).
However, no analysis model exists for facial growth, and the classification of facial growth as a visual expression clustering model is limiting. In this study, we analyzed by applying the -assoc option provided by PLINK software, and the results of this analysis are based on statistical models called likelihood ratio test and Wald test. The reason for applying this analysis is that the phenotype we are targeting is not a general quantitative phenotype, but multinomial variables called facial growth pattern. The currently available method for genome-wide analysis of these variables and multiple SNPs was the statistical model provided by PLINK software. Therefore, the significance between the SNP and the phenotype discovered in this study can be understood as an analysis result of whether the SNP has the explanatory power to explain the phenotype. Some genetic studies based on multinomial variables, and among them, we can check an example of applying the same likelihood ratio test as ours 21,22 .
In addition, the number of samples cannot be considered representative of all Koreans. However, this study represents the first attempt to classify the pattern of facial growth, and when data from two independent groups collected at the same time are analyzed and compared, the common result (the frequency of the pattern is more than 70% coincident) overcomes these limitations.
Most facial changes occur before age 18, but growth and facial remodeling have been shown to continue throughout life. The facial skeleton is generally believed to expand continuously throughout life 23 , which is reflected in the gradual increase in certain facial anthropometric measurements with age, such as anterior nasal cavity and facial width. Certain measurements increase significantly with aging, but some measurements are reduced. The chin length becomes shorter as the mandible of the face grows backward due to aging, resulting in a shorter overall face length.
Some extrinsic variables like gender 24 , body mass 25 are known to effect facial morphology. The main influence of gender on facial phenotypes was reported as nasal area and upper facial area, and body mass index (BMI) was reported as a face width characteristic 24 . Obesity-related sites such as cheeks and neck were excluded from the measurement. So, it is thought that the effect of the degree of obesity in this study is relatively small.
GWAS results provide a hypothesis-free approach to identifying important genetic variations that underlie craniofacial shape differences within populations 26 . A total of 97 significant or suggestive SNPs in 19 gene regions and loci that have previously been associated with facial morphology were identified in this study. For 19 loci showing significant and suggestive phenotypic associations, substantial literature was identified associating these loci with facial development, as shown in Fig. 3. In the current work, we found 10 suggestive SNPs in the horizontal region: FOXK1 27 , IGSF10 28 , FAM161A 29 , POU3F2 30 , DYNC111 31 , SFSWAP 32 , TRIM29 33 , RAPGEF1 34 , PCDH7 35 , and CXCR4 36 . We also found 9 suggestive SNPs in the vertical region: ZSWIM6 37 , CSN3 38 , ATXN1 39 , COL18A1 40 , CHST9 41 , CTNNA3 42 , ASTN2 43 , TUSC3 44 , and MTCL1 45 . The gene annotations from the UCSC database (https:// genome. ucsc. edu) was used to predict the functional effects of the variants. The genes reported to affect embryonic development from UCSC database included CXCR4 46 in the horizontal region and CSN3 38 and TUSC3 44 in the vertical region. The genes reported to affect cranial growth and brain development included ZSWIM6, ATXN1 38 , ASTN2 43 , and MTCL1 45 in the vertical region. In addition, genes related to molecular mechanisms in the regulation of skeletal muscle and cartilage included FOXK1 27 , RAPGEF1 34 , and IGSF10 28 in the horizontal region, and COL18A1 40 and CHST9 41 in the vertical region. The genes associated with retinal circuit components and the growth of sensory organs were FAM161A 29 , POU3F2 30 , and SFSWAP 32 in the horizontal region. Finally, the genes related to frontonasal and dysmorphic facial features were PCDH7 35 , TRIM29 33 , and DYNC111 31 in the horizontal region and CTNNA3 42 in the vertical region.
The authors have generated an interesting report using 2 dimensional images to perform a facial shape GWAS study focused not on fixed time points, but on ontogenetic growth trajectories. To my knowledge, this is a novel analysis. It is also performed in an interesting way, using a mixture of photo types. The argument for this paper is framed around using this for the purpose of being able to improve the accuracy of age progression for missing children, but there is very little discussion of that and this is a very intriguing basic science question.
The facial growth patterns that occur in childhood remain poorly understood and estimating facial growth simply through photographic analysis or the use of existing facial indicators can be difficult. In addition, an individual's unique facial morphologies can be difficult to quantify using simple photographic indicator analysis. The characteristics of facial growth should consider differences in each individual's innate genetic makeup. This study is meaningful because we have classified and characterized facial growth patterns that occur in childhood and will contribute to research on face growth, face recognition, and potentially contribute to finding missing children in the future. This study can serve as a basis for understanding facial morphology and can be expanded to various research fields exploring facial growth, including forensic sciences for both adults and children.

Materials and methods
Study participants. The facial images were obtained from the Human ICT, a company specializing in collecting face data. Related individuals among the participants were not included in the analysis. The facial measurement data were obtained from the National Project of the Missing Child conducted by the Korea Institute of Science and Technology (KIST). Two independent populations (POP1 and POP2) between the ages of 18 and 20 were recruited in different periods: 172 individuals in POP1 were recruited from January 2019 to July 2020, www.nature.com/scientificreports/ and 100 individuals in POP2 were recruited from July 2020 to September 2020. Relatives were not included, only individuals were included. Two types of facial photos were collected: a current picture and older pictures of the same individual. Each individual's current photo was taken in a studio using a Canon EOS 1300D camera with 2592 × 1728 resolution at 400 lx illumination. The participants were asked to close their mouths, hold their faces in a neutral expression, and prevent hair from covering their foreheads. The older photos comprised various ages that can constitute individual chronologies, and a minimum of 4 to a maximum of 21 photos were collected from each participant. Among older facial photos, high-quality photos, such as passport photos and school graduation photos, were prioritized. We also asked the participants to supply photos in which the face was in a neutral expression, not family photos, and torn photos or those that were creased over the face were excluded. Age information was collected for each of the past photos. All facial photographs were digitized with a scanner. We collected a total of 172 current photos and 884 past photos for POP1 and 100 current photos and 600 past photos for POP2. Craniofacial measurements. To extract features for craniofacial measurements in current and past facial photographs, we first detected the facial region in each photograph using Dlib 47 . Dlib 47 detects the face and automatically identifies facial landmarks after the facial region is detected. Among the numerous facial feature points detected, we selected 19 feature points and utilized two facial landmark detectors to automatically extract those feature points, as shown in Fig. 1. One of the detectors was used to extract the facial landmarks using an hourglass network-based feature adaptation network (FAN) 48 approach, whereas the other detector was an inhouse program that combined Dlib 47 and Stasm 49 . The FAN method stacked three hourglass networks, including residual architecture, which is parallel, hierarchical, and multi-scale blocks 50 , to enhance the performance of the feature localization. Stasm 50 is an active shape model (ASM) 51 method with feature descriptors that we fused to the Dlib 47 program. All detectors were programmed in C++ and python in a Qt environment. Facial detection and feature point extraction were performed as automatic processes but could also be manually modified to obtain more accurate feature locations by a well-trained operator (accuracy: 98.8% on average). Euclidean distances between two selected points based on the 19 selected features were calculated, as shown in Supplementary Table 2. A total of 21 facial metric values were calculated. Before this process, the distances www.nature.com/scientificreports/ between the centers of both eyes were normalized to 1 for all images to avoid issues associated with scale differences between components and differences in the Z-axis between the subjects and the camera. Distance calculation programs were implemented in Visual Studio C++.
Face measurement quality controls and sample filtering. We clustered measurements into groups of horizontal and vertical measurements and selected 10 phenotypes to represent the horizontal index and 11 phenotypes to represent the vertical index. The 21 facial phenotypes were measured in each current and past profile photograph, using the 19 facial landmarks in Supplementary Table 2. For all facial phenotypes, we performed data quality controls on volunteers, regardless of the trait being studied. We drew boxplot diagrams for each of the 21 measurements to exclude outlier data from the top 2% and the bottom 2% from all past and current measurements using the R programming language used in R packages 52 . Obesity-related sites such as cheeks and neck were excluded from the measurement. An ordinal multinomial model was applied to show that the results were not due to bias. In this study, we analyzed by applying the -assoc option provided by PLINK software, and the results of this analysis are based on statistical models called likelihood ratio test and Wald test 21,22 .
Early age facial growth pattern (EAFG) analysis. Facial growth shows large variations among individuals; therefore, we graphed the time series of individual facial measurements based on age relative to the current age ( Supplementary Fig. S1) using a non-linear model in QQ plot in R package 53 . We defined a total of 5 EAFG growth patterns that were clustered as follows (Fig. 2): Pattern 1 (DD), continued decrease; Pattern 2 (DI), decrease to increase; Pattern 3 (CC), constant; Pattern 4 (ID), increase to decrease, and Pattern 5 (II), continue increased. Among 21 phenotypes, we coded the 5 EAFG patterns and summarized each individual with similar aging trends in Table 1.
Genotype data. Oral swab samples (KIST) were obtained, and DNA was extracted using ExgeneTM Tissue SV (GeneAll, Seoul, Korea). All DNA samples were amplified and randomly portioned into 25-125-bp fragments, which were purified, resuspended, and hybridized to an Axiom array (TPMRA chip, Thermo fisher, Seoul, Korea), a customized array based on the Asian Precision Medicine Research Array (Thermo fisher Scientific, Waltham, Massachusetts, USA), Following hybridization, the bound targets were washed under stringent conditions to remove non-specific background to minimize noise resulting from random ligation events. The SNP set was filtered based on genotype call rates (≥ 0.98) and MAF (≥ 0.10). Hardy-Weinberg equilibrium (HWE) was calculated for individual SNPs using an exact test. All SNPs reported in this manuscript demonstrated HWE p-values > 0.0001. After filtering, 560,795 polymorphic SNPs were analyzed on chromosomes 1-22.

Imputation of SNPs.
We conducted an imputation analysis to increase the genome coverage. Imputation of genotypes was performed using minimac4 54 at the Michigan Imputation Server (MIS) using the 1000G Phase 3 v4 reference panel 20 . We uploaded phased GWAS genotypes and received imputed genomes in return. After imputation, 7,375,270 polymorphic SNPs were analyzed on chromosomes 1-22. INFO score is over than 0.8.

Genome-wide association studies of EAFG patterns.
To identify not only individual indicators but also indicators that commonly affect facial growth, we performed an analysis that combined the phenotypes of POP1 and POP2. We also conducted a genome-wide association scan of the coded 1 to 5 EAFG growth patterns using asymptotic analyses (likelihood ratio test and Wald test) using the combined population of POP1 and POP2. Population-specific and combined population analyses were performed using PLINK version 1.9 (https:// www. cog-genom ics. org/ plink/) 55 , SPSS (IBM SPSS Statistics Inc., New York, U.S.) 56 , and R Statistical Software 52 . We calculated the beta coefficient and the standard error (SE) values for the association study. To compare the GWAS results for each population, we conducted a replication study using 172 samples from POP1 and 100 samples from POP2. We selected the genetic markers associated with the 5 EAFG patterns in each GWAS, determined by association p-values < 1 × 10 −5 in the combined dataset and p-values < 0.05 for the individual population datasets and the replication study. In this study, GWAS analysis was performed based on about 800,000 SNPs. The Bonferroni correction p-value threshold is applied. The results are shown in Tables 2 and 3 and Manhattan plots are depicted in Supplementary Fig. S3. The QQ plot, generated using R Statistical Software 52 , of the observed p-values showed minimal inflation of the GWAS results from the combined population sample ( Supplementary Fig. S2).

Annotation of SNP-associated genes.
To identify and annotate genes that are functionally related to suggestive and significant SNPs identified in the GWAS, SNP locus data were obtained from the UCSC Genome Browser (Genome Bioinformatics Group, University of Santa Cruz, Santa Cruz, CA, USA). The gene annotations from the UCSC database and Genotype-Tissue Expression (GTEx) database (GTEx Analysis Release v.8, http:// www. gtexp ortal. org/) were used to predict the functional effects of the variants. We constructed regional plots of association for regions of interest using the program LocusZoom ( Supplementary Fig. S4) 57 .

Data availability
Raw genotype or phenotypic data cannot be used due to limitations imposed by ethics. The Summary statistics obtained here are based on the GWAS analysis and can be accessed with the supplementary materials.